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ABSTRACT 

We describe results from time-dependent numerical modeling of the collisionless reverse shock terminating 
the pulsar wind in the Crab Nebula. We treat the upstream relativistic wind as composed of ions and electron- 
positron plasma embedded in a toroidal magnetic field, flowing radially outward from the pulsar in a sector 
around the rotational equator The relativistic cyclotron instability of the ion gyrational orbit downstream of 
the leading shock in the electron-positron pairs launches outward propagating magnetosonic waves. Because 
of the fresh supply of ions crossing the shock, this time-dependent process achieves a limit-cycle, in which the 
waves are launched with periodicity on the order of the ion Larmor time. Compressions in the magnetic field 
and pair density associated with these waves, as well as their propagation speed, semi-quantitatively reproduce 
the behavior of the wisp and ring features described in recent observations obtained using the Hubble Space 
Telescope and the Chandra X-Ray Observatory. By selecting the parameters of the ion orbits to fit the spatial 
separation of the wisps, we predict the period of time variability of the wisps that is consistent with the data. 
When coupled with a mechanism for non-thermal acceleration of the pairs, the compressions in the magnetic 
field and plasma density associated with the optical wisp structure naturally account for the location of X-ray 
features in the Crab. We also discuss the origin of the high energy ions and their acceleration in the equatorial 
current sheet of the pulsar wind. 

Subject headings: acceleration of particles - ISM:individual (Crab Nebula) - pulsars:general - pulsars: indi- 
vidual (Crab Pulsar) - shock waves 



1. INTRODUCTION 

The mechanisms through which compact objects excite sur- 
rounding synchrotron nebulae are one of the long-standing 
problems of high-energy and relativistic astrophysics. The 
loss of rotational energy from the central magnetized ob- 
ject underlies the synchrotron emission observed from pulsar 
wind nebulae (PWNs; e.g., Bandiera, Amato, & Woltjer 1998) 
and is one of the prime candidates for the energization of the 
jets in active galactic nuclei, radio galaxies (see Krolik 1999), 
and gamma ray bursts (e.g., Blandford 2002). PWNs form the 
nearest at hand laboratories for the investigation of the physics 
of such high-energy particle acceleration, and are the example 
where rotational energization clearly bears responsibility for 
the nonthermal activity. 

Despite intensive study, however, the physics of how the 
rotational energy gets converted to the observed synchrotron 
emission has not been satisfactorily understood. Magnetohy- 
drodynamic wind models (Kennel and Coroniti 1984a,b and 
many subsequent efforts) provide an adequate macroscopic 
description of the "bubbles" of synchrotron-emitting particles 
and fields surrounding young (t ^ Kfi yr) rotation-powered 
pulsars, but identification and quantitative modeling of the 
mechanisms that convert the unseen energy flowing out from 
the central pulsar into the visible synchrotron emission have 
remained elusive. 

Most attention has focused on the idea that the outflow is a 
supersonic and super-Alfvenic wind, perhaps with embedded 
large-amplitude electromagnetic waves (Rees and Gunn 1974; 
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Melatos and Melrose 1996). Kennel and Coroniti's model 
simplified the picture by assuming the wind to be steady, with 
no wave structure. In that case, most models attribute the con- 
version of relativistic outflow energy to the nonthermal parti- 
cle spectra that emit the observed synchrotron radiation in the 
body of the PWN to a relativistic collisionless shock wave 
that terminates the flow. This idea was originated by Rees 
and Gunn, who also observed that the shock should form at 
a radius where the dynamic pressure of the outflow is ap- 
proximately equal to the pressure of the surrounding nebula. 
They also noted that in the case of the Crab Nebula, this ra- 
dius approximately corresponds to the location of the moving 
"wisps" discovered 80 years ago (Lampland 1921) in optical 
observations of that PWN. This identification has motivated 
many studies of the wisps, as possible direct manifestations 
of the shock wave thought to underlie the transformation of 
outflow energy into the nonthermally emitting particles in the 
nebula. 

Shortly after the discovery of the Crab pulsar, Scargle 
(1969) described these motions in more detail and attempted 
to interpret them as wave phenomena in the nebula stimulated 
directly by glitches in the pulsar's rotational spin-down (Scar- 
gle & Harlan 1970; Scai-gle & Pacini 1971). Subsequent ob- 
servations made clear that the association of wisp activity to 
pulsar rotational activity was a spurious consequence of un- 
dersampling the variability of the wisp motions. 

The advent of the magnetohydrodynamic model has fo- 
cused attention on the wisps' structure and variability as a 
probe of the physics of the region where the pulsar's outflow 
interacts with the nebular plasma. Substantially improved ob- 
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servations in the optical (Hester et al. 1995; Tanvir, Thomson, 
& TsikarishviH 1997; Hester 1998a; Hester et al., 2002), ra- 
dio (Bietenholtz, Frail, & Hester 2001) and X-ray (Mori et 
al. 2002; Hester et al. 2002) now make it possible to test 
relatively refined theories of the shock's structure and particle 
acceleration properties. 

A number of ideas have been suggested for the physics be- 
hind the observed dynamic behavior. Hester (1998a) sug- 
gested that the wisp structures reflect synchrotron cooling 
at constant pressure of the particles whose outflow mo- 
menta randomize at a shock in a magnetized electron-positron 
plasma to a power-law distribution of particles with slope in 
energy space flatter than E"^, an idea that requires the X-ray 
spectrum of the plasma in the wisp region to be flatter than 
e~''^ (X-ray flux in photons/keV-cm^-sec). However, only the 
particles with energy high enough to emit nebular gamma rays 
(e > lOMeV) have synchrotron loss times short enough to 
create a cooling instability whose time scales are compara- 
ble to the observed variations in the optical and X-ray, and 
the gamma ray spectrum of the Crab (presumably arising in 
a compact region around the pulsar no larger than the X-ray 
source) has a spectral index steeper than 1.5. Chedia et al. 
(1997) suggest that a drift instability of a subsonically expand- 
ing plasma (no shock wave) can explain the wisp structure - 
how the outflow from the pulsar escapes catastrophic adia- 
batic losses is not addressed, however Begelman (1999) sug- 
gests that the wisps represent the nonlinear phase of a Kelvin- 
Helmholtz instability between subsonic layers expanding in 
the rotational equator of the pulsar and flow at higher abso- 
lute latitude moving with a different velocity. 

All of these ideas follow Kennel and Coroniti in assuming 
that the pulsar's outflow decelerates in a shock wave of un- 
observably small thickness, with the wisp dynamics due to 
phenomena in the post shock plasma (or, in the case of the 
Chedia et al. model, do without a shock completely). In con- 
trast. Gallant and Arons (1994, hereafter called GA), building 
on earlier work on the structure and particle acceleration char- 
acteristics of relativistic shock waves in collisionless, magne- 
tized electron-positron (Langdon, Arons, & Maxl988; Gal- 
lant et al. 1992) and electron-positron-ion (Hoshino et al., 
1992) plasmas, attributed the observed wisp structure to the 
internal structure of an electron-positron-ion shock wave. The 
idea is motivated by the fact that the gyration of the heavy ions 
within the shock transfers energy to positrons and electrons 
through emission of ion cyclotron waves (formally, the mag- 
netosonic mode of the magnetized electron-positron plasma). 
This leads to the formation of nonthermal particle spectra in 
rough agreement with the particle spectra required to explain 
the optical. X-ray, and gamma-ray emission, if the ions carry 
a large fraction of the pulsar's total energy loss, at least in the 
equatorial sector illustrated in Figure^ The deflection of the 
ions' outflow by the magnetic field deposits a large amount 
of compressional momentum in the magnetic field, which is 
transverse to the flow, and into the plasma frozen to the 
field. These compressions appear as brightenings in the sur- 
face brightness of the nebula. GA showed that their steady 
flow model gives a good account of a single high-resolution I- 
band image of the wisps (van den Bergh & Pritchett 1989). 
However, the wisps are known to be highly variable, on a 
time-scale of weeks to months. An adequate account of their 
behavior requires a time-dependent theory. 

In this paper we present a model for time-variability of the 
wisp region in the Crab by extending the work of GA to in- 
clude the time dependence in the flow. As was realized early 



on by Hoshino et al. (1992) and Gallant et al. (1992), the 
coherent ion orbits invoked in the model of GA are unstable 
to gyrophase bunching and emission of magnetosonic waves. 
We find that this instability introduces a time dependence into 
the shock in the pair-ion plasma that causes the variable sur- 
face brightness to have substantial similarity to the observed 
variability of the wisps. Because of the presence of the con- 
tinuous stream of fresh ions passing through the pair shock, 
the relativistic ion cyclotron instability provides a mechanism 
for sustained periodicity and wave emission dynamics that 
closely reproduces time-resolved observations of the wisps. 

In ^we describe the model and the underlying assump- 
tions. In §^ and 13 we present the results of the simula- 
tions and the dynamical predictions for the wisp region and 
compare them to observations. In ^ we discuss the model's 
successes, with emphasis on its applicability to other PWNs 
that have wind-nebula dynamics not very affected by radiation 
losses, as well as the model's limitations, pointing out fur- 
ther improvements needed in order to make fully quantitative 
comparisons with observations. We also discuss some obser- 
vational tests of th e ide as at their current level of development 
in this section. In ^5.4l we describe a scenario for the origin of 
the high-energy ions and their acceleration in the equatorial 
region of the pulsar wind. Our conclusions are summarized 
in ^ The technical details of our hybrid numerical approach 
are left for the Appendix. 

2. THE MODEL 

Following GA, we assume that the equatorial wind of the 
Crab pulsar consists of relativistic electron-positron pairs and 
an ionic component, which is minor in number density but 
energetically significant. Because the ion gyration timescale 
is much larger than that of the pairs, we resort to a hybrid 
model of the termination shock in the wind. In it, the ions are 
treated kinetically, while the pair component is modeled as 
an MHD fluid. When the magnetization of the wind is small, 
the system is nearly charge and current neutral (see Appendix 
for derivation; eq. JA14IA15t ). Therefore, the charge and 
current density in the pairs can be found by following the ion 
dynamics. 

We assume the flow to be spherically symmetric in a 20° 
sector about the equatorial plane. Using the spherical coordi- 
nates (r, 6, (p) centered on the pulsar, and allowing for vari- 
ation of all quantities with r only, the conservation laws for 
mass, momentum, energy, and the field evolution equations 
can be written as 

|(P7)+;^|:('-W.) = 0, (1) 

at or or r 

--phj^(3l + eNiiPw - (3g)B^ = 0, (2) 
r 

^(phj^Pe) + \^(r^phj^f3rPe) + -(phj^Me) 
Ot r^ or r 

+emj3r-l3ir)B^ = 0, (3) 

O \ d 

-(phj^-Pi_)+-—(r^phj^Pr) 
ot r^ Or 

+em(3ief3r-f3irf39)B4, = Q, (4) 
d 19, B^Br 

-B^ + -—(r^B^Pr) ^ = (5) 

ot r^ or r 

where p and h are the pair density and specific enthalpy in 
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Fig. 1 . — (a) HST optical image of the inner Crab Nebula (after Hester 1995), (b) Zoom in on the wisp region, (c) X-ray nebula from Chandra (after Weisskopf 
2000; [b] and [c] are nearly to scale); (d) Illustration of the polar jet and the equatorial flow geometry of the Crab wind. This figure is also available in color in 
the electronic edition of the Astrophysical Journal 



the proper frame, 7/3 is the 4-velocity of the pairs, A^, = n,7, 
is the observer's frame density of the ions, and is the 
toroidal component of magnetic field, which is the only mag- 
netic field component that we consider in the equatorial plane. 
P±^ and P|| are the components of the pair plasma's pressure 
tensor perpendicular and parallel to B. The field is normal- 
ized by the upstream value B\, and expression (|6j defines e. 
Since we are interested in the dynamics on the ion gyration 
scale, we normalize time by Lo'j^ and radius by ru\, where 
u)ci\ = ZeBi/'ynniiC and rm = c/ucn are the upstream values 
of the relativistic cyclotron frequency and ion Larmor radius, 
respectively, and Ze is the ion charge. We normalize the en- 
thalpy and pressure by the upstream kinetic energy of the pairs 
m±c^A^i±7i, where A^i± =^in\± is the pair number density in 
the observer's frame. After traversing the shock, the pairs are 
confined to the plane perpendicular to the toroidal magnetic 
field until pitch-angle scattering isotropizes the pair distribu- 
tion function. We allow for the anisotropic pressure in the 
directions perpendicular and parallel to the magnetic field, 
P±^ P\\- These components of the pressure tensor are re- 
lated through the adiabatic index F such that P± = (F- 1)£, 
where e = + P|| is the internal thermal energy of the pairs. 
This results in F = (3 + P||/P_i)/(2 + P||/P_i) with F decreas- 
ing from 3/2 to 4/3 as pressure isotropizes. In our model 
we simulate the effect of pitch-angle scattering in the pairs 
by gradually decreasing F with distance from the pair shock 
over a distance of 0.05 pc (for further discussion of theory and 
simulation of pitch-angle scattering in magnetized electron- 
positron-ion plasmas see Yang et al. 1993; Yang, Arons, & 
Langdon 1994). The relative fractions of the preshock wind 
energy carried by the ions, the pairs and the Poynting flux are 
parameterized using 
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where a± and cr,- are the initial magnetization of the pairs 
and ions, respectively. Global models of the nebula (e.g.. 
Kennel and Coroniti 1984a) constrain the total magnetization 
(Ttot = (l/ci + l/cr,)~' to be small, ~ 0.003. Because of the 
small initial magnetization of the wind and our attention to 
the region near the termination shock, we neglect the terms 
proportional to cr that otherwise would enter into ( I1I4> (see 
Appendix for the full equations). 

Equations (Ill4t describe the system of evolution equations 
that we solve numerically to find the time- variability of the 
colUsionless shock structure. We use the particle-in-cell (PIC) 



method (Birdsall and Langdon 1991) for advancing ions and 
calculating their current and charge density, while the pairs 
are evolved by finite-differencing the MHD equations on a 
one-dimensional mesh. Particular attention is paid to the 
boundary conditions, with an outflow nonreflecting condition 
(Thompson 1987) implemented on the downstream boundary. 

The inner boundary condition for the pair plasma and mag- 
netic field comes from the jump conditions for an infinitesi- 
mally thin shock wave in the pair plasma, as in Kennel and 
Coroniti (1984a). The shock is located inside the domain and 
the flow is extrapolated upstream to the injection boundary at 
0. 1 ru. The ions enter as a cold stream at the inner boundary, 
all with identical energy 71 m,c^, where 71 is the Lorentz fac- 
tor of the upstream wind - we assume that the upstream ions 
and pairs flow with the same 4-velocity. The jump conditions 
at the pair shock are such that the magnetic field increases by 
a factor of 1/(F- 1) over its value in the upstream wind, and 
the velocity of pair flow decreases by the same factor, while 
the tangential electric field is continuous through the shock. 
The ion stream, therefore, starts to gyrate as well as E x B 
drift in the magnetized pairs as shown in Fig. 12^. The turning 
points in the ion orbits are regions of strong vertical current 
that, through the j x B force, provide a compression on the 
pairs. This compression results in locally higher density of 
the pairs at each ion turning point, and compression of the 
frozen-in magnetic fields. This field, in turn, affects the mo- 
tion of the ions, and hence the initial configuration is bound to 
change. Although it is possible to load the ions and the fields 
in a self-consistent static configuration, we find that the result 
of the dynamical simulation is independent of the initial load 
method. We describe the dynamical simulations in the next 
section. 

3. SHOCK DYNAMICS 

The coherent ion loops generated during ion transit through 
the pair shock are an excellent environment for the growth of 
a relativistic ion cyclotron instability. This is a variant of the 
classic instability of a "ring in momentum space" and pro- 
ceeds via gyrophase bunching and subsequent phase mixing. 
The crucial distinction in the shock case from the develop- 
ment of the instability in a uniform plasma is that there is a 
preferred direction of the flow, so that the time development of 
the instability is spread out in space as the ions undergo E x B 
drift into the shocked pairs. The typical growth rate is on the 
order of the ion Larmor time (Hoshino & Arons 1991) in the 
postshock field, and within the first ion loop gyrobunching 
develops characteristic "knots," shown in FigurelJl^. Beyond 
several orbits in the postshock field the instability has taken its 
course, and the downstream ion orbits become chaotic, with 
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0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 0.5 1 1.5 2 
r/ru, r/rLii r/r^i, r/fu, r/ri,, r/r^, 

Fig. 2. — Phases of the limit cycle. Top: Magnetic field, in units of preshock value B\, vs. distance. Middle: Ions' momenta through the shock, with ion 
momentum in units of the upstream momentum vs. distance, with length units of preshock ion Larmor radii. Bottom: Ion density in the radial direction, with 
density normalized to the preshock ion density. The unit of length is the ion Larmor radius based on upstream parameters, rm. The columns are as follows: 
(a) Initial state of the F = 3/2 shock with F decreasing to 4/3 by r = \.5rm . The cold ion stream is loaded through the pair shock without affecting the field 
structure, (b-f) Developed limit cycle. A gyroknot advects around the ion ring and triggers the next cycle upon leaving the region. This figure is also available as 
an mpeg animation in the electronic edition of the Astrophysical Journal. 



the ion distribution eventually relaxing to a quasi-Maxwellian 
form. However, in the case of a shock in a relativistic flow 
and in our simulation, the cold ions are continuously flowing 
through the shock, thus constantly refreshing the instability's 
flagging spirits. 

The interplay between the low entropy input of ions and the 
randomizing cyclotron instability creates interesting dynami- 
cal consequences for the region of the first ion gyration loop. 
In this region we observe the development of a limit cycle, 
in which the first ion loop develops gyrobunched ion knots 
that are advected downstream with the overall E x B flow. As 
one knot is leaving this region, it creates an enhanced mag- 
netic field that triggers the development of a new "gyroknot," 
which goes around the ion ring and the process repeats. The 
snapshots of the phases in this limit cycle are shown in Fig. 

In these frames we display the magnetic field, the ion 
stream flowing through the shock, and the corresponding ion 
density. The postshock magnetic field shows the character- 
istic "double-horn" structure due to compressions in the pair 
density at the turning points of the zeroth-order ion drift or- 
bit. The development of the cyclotron instability in the first 
ion loop introduces additional compression in the field be- 
cause the gyrobunches in the ion stream exert additional j x B 
stress on the pairs. As the gyroknot is advected around the ion 
loop, the associated compression in the pairs launches magne- 
tosonic waves. The area of the first ion ring can be thought of 
as the near zone of a radiating region, while the flow further 
downstream represents the wave zone. As the wave propa- 



gates between the turning points in the first ion loop, it im- 
prints a moving component in the double-horn structure of 
the magnetic field. When this compression in the field reaches 
the second ion turning point, it triggers gyrobunching in the 
incoming ions. The period of the cycle is found from simu- 
lations to be ^ 0.5Tu, where Tu = Itt/lou is the Larmor time 
for the ions in the postshock field. The cycle period can be 
understood as roughly the time in the laboratory frame that 
it takes for the ion to go around one loop in the E x B drift 
orbit, or [(1 -/3j)/(l + Pd)]^^^Tu, where (3dC = (T-l)c is the 
initial drift speed in the downstream region. 

4. APPLICATION TO THE WISPS' DYNAMICS 

As did GA, we interpret the wisp region in the Crab as a col- 
lisionless shock in the pair plus ion plasma. Although the ions 
are not directly observed, the compressions in the pair den- 
sity and magnetic field associated with the turning points in 
the ion orbit should result in regions of increased synchrotron 
brightness. GA associated the first two ion turning points with 
wisp 1 and wisp 2 (see Fig. ^). Our time-dependent analysis 
augments this interpretation with additional features. We find 
that the back-and-forth motion of the gyro-knot in the region 
of the first ion ring acts as a periodic source of magnetosonic 
waves that make the shock structure very dynamic. The mag- 
netosonic waves that propagate back toward the pair shock 
are reflected from it and compress the magnetic field near the 
shock, brightening wisp 1 . As the waves propagate from wisp 
1 to wisp 2, they appear as a moving bright feature, a "moving 
wisp," that brightens the region of wisp 2 when it crosses it. 
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Fig. 3. — Top: //5r snapshots of the wisp region in the Crab. Separation between the first and the last pair of frames is 3 weeks; middle frames ai'e separated by 
a week. The observation interval is from 1995 December 29 to 1996 February 22. Data were extracted from the HST archive (observations by J. Hester). Bottom: 
Simulated suii'ace brightness of the wisp region. The moving wisp is being emitted from the first wisp and propagates toward the second wisp. Fine filamentary 
structure is created in the second wisp region. This figure is also available in color and as mpeg animation in the electronic edition of the Astrophysical Journal. 



The sequence of propagation and brightening of the features is 
choreographed by the underlying Hmit cycle of the instability 
in the first ion loop. The overall timescale of the limit cycle 
can be obtained by choosing the preshock ion Larmor radius 
ru\ such as to fit the observed separation of the wisps 1 and 
2 with the separation of compressions in the model. The limit 
cycle then occurs with the period Tcyde ~ \Tu = ■K(T—\)rui/c. 
Here we used the magnetic field jump conditions to express 
the postshock Larmor time in terms of preshock quantities. 
For typical values F = 3/2, rm =0.1 pc we get Tcyde ~ 6 
months. 

We have taken our one-dimensional profiles of the magnetic 
fields, pair density and pair velocity to construct the sequence 
of simulated surface brightness maps of the wisp region in the 
optical band. We assume that in the proper frame the pair 
flow emits synchrotron radiation as a relativistic Maxwellian. 
This is essentially a monoenergetic emission and represents 
the surface brightness distribution in space for any relatively 
narrow bandpass from pairs whose Larmor radii are small 
compared to the ion Larmor radius (photon energies below 10 
MeV). The radiation from the flow moving along the equato- 
rial plane, tilted by ^ 65° from the plane of the sky, is further 
Doppler boosted and beamed on the way to the observer The 
resulting maps, shown in Fig. |3l are a qualitative picture of 
the way the features in the wisp region should appear The 
Doppler beaming makes the faster flow near the shock ap- 
pear to have smaller angular extent than slower flow in the re- 
gion of wisp 2. For comparison, we also reproduce snapshots 
from 1995-96 Hubble Space Telescope (HST) observations of 
the inner nebula during a moving feature emission episode 
in the same figure (Hester et al. 1998b). Observations from 
the 2000-2001 campaign (Hester et al. 2002) show very 
similar dynamics of the features. In principle, the simulated 
brightness profile can be constructed for the X-ray band as 
well and will look very similar to the optical image. This pro- 
cedure will, however, overpredict the amount of beaming and 
make the inner Chandra ring appear as two disjoint arcs, con- 
trary to observations. We comment further on this discrepancy 



in 1521 

Using the simulated surface brightness maps we can de- 
scribe the time dependent behavior of the wisp region 
throughout the cycle of emission according to our model: 

1) The emission of a major moving wisp is preceded by 
a brightening and thickening of wisp 1 . This may also 
appear as a forward motion of wisp 1 . After the emis- 
sion, it dims and retracts to the original position. In 
general, the oscillating pressure of magnetosonic waves 
in the wisp region causes oscillation in the position of 
the shock in the pairs. This can make wisp 1 appear as 
oscillating with the cycles of emission of the moving 
wisps. The amplitude of this oscillation is relatively 
small (< 10% of the distance between wisp 1 and the 
pulsar). 

2) The compression associated with the moving wisp in- 
creases in brightness while propagating toward wisp 2. 
Because of this, the moving wisp appears to turn on 
some distance away from wisp 1 (< 20% of the dis- 
tance between wisp 1 and wisp 2). The moving wisp 
appears as a fairly sharp feature. 

3) The region of wisp 2 is located at the interface of the 
first and second ion loops. By the second gyration loop 
the ion cyclotron instability is well into the nonlinear 
stage, and the ion orbits are scrambled and chaotic. 
This ion background provides a large amount of short- 
wavelength features that comove with the flow as the 
ions drift downstream. In fact, observations by Hester 
et al. (1995) indicate that wisp 2 has a fine fibrous struc- 
ture. We attribute the creation of such filaments to the 
chaotic ion background. Although the fine filaments 
should appear moving away from the shock roughly 
with the flow velocity, we expect that new filaments can 
appear in this region without having a precursor travel- 
ing all the way from the shock. 
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4) When the moving wisp propagates toward wisp 2, it 
brightens up wisp 2 and its surroundings. The mov- 
ing wisp is a compression in the field and pair density 
that propagates as a fast magnetosonic wave. This wave 
was emitted from the first postshock ion gyration loop 
and travels at the fast magnetosonic speed on top of the 
background pair flow. In the case of low magnetiza- 
tion the wave speed is given by c(/3flow+ Vr- 1)/(1 + 
/3flowVr- 1). Assuming that the pairs have isotropized 
in the region of the second wisp, so that F = 4/3, we get 
Vwisp ~ 0.6c as would be seen in the equatorial plane of 
the Crab. When the moving wisp traverses wisp 2, it 
overtakes slower moving fibrous filaments, and causes 
more filaments to be created in its wake, as it scatters 
on the random ion background. 

5) The propagation of the fast-moving wisp from the first 
to the second wisp takes about half the cycle of the wisp 
emission 90 days). For the second half of the cycle 
no major brightening occurs, until the first wisp starts 
to brighten again, indicating the beginning of the next 
cycle. This behavior paradoxically results in the wisp 
structure appearing both stationary and moving at the 
same time. On average there are always two bright re- 
gions corresponding to the first and second wisps, with 
the second wisp being more diffuse. These regions cor- 
respond to compressions due to the turning points in the 
first ion loop, which are more or less fixed in space, but 
can oscillate in brightness. At the same time, the in- 
stability in the ion gyration produces intermittent wisps 
that move between the first and second wisp at the mag- 
netosonic speed and change the brightness of the main 
wisps as they compress the plasma and the field. Yet, 
the first and second wisps do not move as a wave them- 
selves. 

The model has several parameters that we can vary to study 
the shock dynamics in qualitatively different regimes. The 
main parameters are the ratio / of the preshock ion Larmor 
radius to the radius of the shock in the pairs, / = rui/r,, 
and the relative fraction e of the energy density carried by the 
ions and the pairs (eq. (6l). The parameter / controls the 
degree to which spherical geometry affects the shock dynam- 
ics. For small values of /, the background flow is essentially 
plane-parallel on the ion gyrational scale, and the ions do not 
experience the gradient in the background magnetic field. For 
the Crab parameters, f ^ 1, so the effects of spherical ge- 
ometry are important. Increasing / decreases the separation 
of the turning points in the ion orbit and changes the relative 
strength of plasma compressions corresponding to the turning 
points, as illustrated in Figure 7 of GA. The parameter e con- 
trols the amount of influence the ions have on the pair fluid. 
For small e the underlying pair shock structure is hardly per- 
turbed. This makes the instability grow at a slower rate, and 
the first several ion loops may appear steady. The gyro-knots, 
characteristic of the instability, then appear only several Lar- 
mor radii downstream. In the limit e there is, of course, 
no instability, and the ions and the pairs are decoupled. For 
finite e, once the instability proceeds, the period of the limit 
cycle is not sensitive to e, and is approximately half of the 
local Larmor time of the ions. We fix the value of / « L2 to 
match the separation of the wisps with the compressions at the 
turning points in the ion orbits. We then vary e to achieve the 
conditions where both the instability grows sufficiently by the 



time the ions cross the first ion loop, and the compressions in 
the magnetic field and the density result in the intensity vari- 
ations consistent with observed variation in the nebula. For 
perturbations 5B/B ^ 1 we need e '--^ 1. For the pair injection 
rate of .;V± « lO^*** s"' needed to explain the X-ray source in 
the Crab, we then find that the ion injection rate should be 
.^V; = em±/miN± « 10^"^ s"^ if we assume that the ions are he- 
lium (the model favors ions with a Z/A ratio of ^ 2, which 
could be either He or Fe). These ions travel with the same 
Lorentz factor (10*) as the equatorial pair wind. 

5. DISCUSSION 

5.1. Model Successes 

The theory outlined here has a number of virtues in its 
description of the physics of the termination of the Crab 
pulsar's equatorial wind. The timescale and length scales 
of the observed optical and X-ray variations can be semi- 
quantitatively reproduced, for entirely reasonable values of 
the ion energy. Granted the approximate validity of the un- 
derlying Kennel and Coroniti (1984a) flow model, for ions 
that accelerate through a fraction ^ 0.1-0.2 of the available 
open field potential $ w 4 x 10"' V to energies ~ 0.15Ze<l> w 
6 X 10'^ eV, the ion Larmor radius, which sets the wave- 
length scale of the moving features, is ru ~ IQ^^B'^f^irJ r) 
cm. Here < 1.5 x 10'^ cm is the radius of the shock in 
the pairs, located at or interior to the Chandra X-ray ring and 

Bs = <T^^^B„ebiiia ~ 16-S.S16 /^G. The ion energy is obtained by 
asking the model to replicate the spacing of the moving fea- 
tures in the optical and X-ray images. The ions are deflected 
and become magnetically well coupled to the pair plasma in 
the increasing magnetic field at r > r,urnmg = [fu(rs)rs]^^^ ~ 
4 X 10'^(7/6.5/-B.si6)'''^ cm. This is the radius at which the 
compressional limit cycle formed by the first ion loop stands 
as a quasi-stationary feature in the pair flow. Given the close 
resemblance between this feature of the model and the behav- 
ior of the Chandra X-ray ring, it is tempting and natural to 
identify the X-ray ring with the first ion loop's turning points. 
If so, the quasi-stationarity of this feature then has a natural 
explanation, as being the consequence of the ion cyclotron in- 
stability within the shock structure, without having to invoke 
extra physics and an otherwise unexplained external stimulus 
to account for the time scale of the ring's repetitive variations. 
The non-axisymmetric knotted structure of the ring might be 
explainable as the consequence of instabilities of the mirror 
type, driven by the reflected ions' pitch angle anisotropy with 
respect to the toroidal magnetic field - a topic well beyond the 
scope of the current investigation, however. 

The model also gives a natural interpretation of the moving 
features in both optical and X-ray, as compressional waves in 
the pair plasma emitted by the time variable ion currents in 
the first loop limit cycle. These modes are a mixture of mag- 
netosonic waves moving with respect to the plasma flow and 
entropy and slow modes frozen into the flow. This behavior 
supports the identification of the Chandra ring as the location 
of the ion limit cycle structure, since that ring's appearance 
suggests it to be the wave emission zone. 

The theory is adiabatic with respect to radiation losses from 
the pairs; that is, the formation of the ring and the moving 
wisps does not depend on the synchrotron loss times of the 
pairs being comparable to or shorter than the local flow time 
in the vicinity of the ring. Since ringlike and wisp features 
have been discovered near pulsars in other PWNs at radii 
where dynamic pressure balance arguments suggest formation 
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of wind termination shocks, yet all occur in systems whose 
synchrotron cooling times even for X-ray-emitting particles 
are larger than the local flow times, the theory developed here 
appears to have general applicability, in contrast to models 
that depend on the unusually strong magnetic fields and short 
synchrotron cooling times prevalent in the Crab Nebula (e.g., 
Hester 1998a; Hester et al. 2002). Gaensler et al. (2002) give 
an example of our theory's use, in the interpretation of the 
X-ray features observed close to PSR 1509-58. 

This model has a substantial theoretical virtue: the required 
ion flow has the magnitude expected if it is the electric re- 
tum current (Goldreich-Julian current) that should exist in the 
rotational equator of the pulsar's wind (Michel 1974; Con- 
topoulos et al. 1999), which prevents the star from charging 
up when non- vacuum torques are a significant contributor to 
spindown. This conclusion is the same as was reached by GA 
from a time-independent flow model, now seen as a robust 
consequence of the fully time-dependent theory. 

5.2. Model Limitations and Failures 

There is a substantial (~20%) discrepancy between the ve- 
locities of the moving magnetosonic wave features computed 
in the model (0.6c) and the speeds of the observed features 
0.5c, in the geometry adopted by Hester et al. (2002), 
itself a velocity large compared to what one would expect, 
{c/2>){Rshock/r)^, for features frozen into a simple Kennel & 
Coroniti 1984a flow]. This discrepancy may be related to 
the geometrical simplifications of the model, discussed fur- 
ther below, or it may be an intrinsic difficulty. 

The model's one-dimensionality is a serious Umitation. We 
assumed the flow to be strictly radial, within a spherical sec- 
tor around the rotational equator Within that sector, we as- 
sumed the magnetic field to be strictly toroidal and have a 
single sense of winding with respect to the rotation axis of the 
pulsar, independent of rotational latitude. Going along with 
this one-dimensionaUzation, we assumed the flow of pairs and 
ions to be completely charge neutral and to be current neutral 
in the radial direction. Yet we concluded that the ion flux 
corresponds in magnitude to an electric current flowing out 
in the equatorial plane with strength sufficient to qualitatively 
and quantitatively alter the magnetic field. 

A more realistic model would take this result seriously and 
treat the ion flux as a net electric current in the radial direc- 
tion (balanced in a quasi-steady state by a polar electron cur- 
rent feeding the polar jets). Such a current in the upstream 
wind separates two hemispheres of toroidal (and very weak 
radial) magnetic field that are oppositely directed in the two 
hemispheres. In such a field, high-energy ions in the equato- 
rial current flow, upon crossing the shock in the pairs, are not 
fully reflected by the compressed magnetic field in the heated 
pair plasma. They are deflected from their initial radial or- 
bits, but transfer a smaller fraction of their outward momen- 
tum into compressions of the pair plasma. Thus, the model 
described here somewhat overestimates the compression of 
the magnetic field in the region where the limit cycle appears, 
probably by a factor on the order of 2. Generating synthetic 
surface brightness maps, as a function of assumed ion and pair 
energy and number fluxes, that are quantitatively reliable re- 
quires taking this more realistic geometry and ion dynamics 
into account. 

Several comments are in order regarding the speed of prop- 
agation of the wisps. Quasi-periodic temporal drive due to 
the ion instability leads us to interpret the moving features in 
the Crab as the fast MHD waves. This, however, leads not 



only to overestimation of the speed of the features as outlined 
above, but also to an uncomfortable prediction that the wisps 
should reach an asymptotic velocity as they propagate down- 
stream. As the background flow decelerates, the speed of the 
disturbance decreases to the sound speed Vr- Ic = 0.58c for 
r = 4/3. Fast waves are not the only type of perturbations 
generated by the instability, however. Ion motion also excites 
entropy, Alfven and slow waves. All of these perturbations 
have zero phase velocity for the toroidal background mag- 
netic field; therefore, they are advected with the flow. The 
entropy wave does not have any perturbation in the magnetic 
field and therefore will not appear as bright as the modes that 
perturb both the density and the magnetic field. The Alfven 
mode is peculiar because it is the only disturbance that has the 
transverse velocity perturbations. It can be advected with the 
ion flow that generates the perturbation in the pairs' latitudinal 
velocity. 

As the fast wave from the emission region crosses wisp 2 
and encounters the chaotic ion background beyond the first 
ion E X B loop, there will be some compressions in the pairs 
and the field that move with the average E x B velocity of 
the ions (which is the same as the pair flow velocity). These 
compressions are not free-propagating, but driven by j x B 
forces from ion streams. They are similar to the double- 
horn soliton structure generated by ion loops (Alsop & Arons 
1988), but they are moving as the ion loops are advected 
with the unstable ion stream. Therefore, at larger distances 
it is reasonable to expect wisplets that appear to move with 
the background flow. However, such features would be rel- 
atively short lived (approximately a month) and would ap- 
pear to come seemingly from nowhere. Therefore, beyond 
wisp 2 we expect to see two populations of wisps - ones 
moving with the background decelerating flow (slow magne- 
tosonic and Alfven waves), and fast magnetosonic perturba- 
tions that can be traced back to the shock in the pairs at wisp 
1 . The intensity of the fast compressions will get diminished 
with distance, and they will naturally blend into slower back- 
ground flow. A clear observational tracking of a single feature 
through the whole flow (with snapshots separated by no more 
than 2 weeks) would be very beneficial to settling the issue 
of the speed of the flow and the nature of the speed of the 
features. 

The front-to-back asymmetry of the optical wisps' bright- 
ness has an obvious interpretation in the Doppler boosting en- 
forced by the underlying flow, which preserves on average 
the Kennel and Coroniti radial profile v,. = (T - \)c(r/Rs)^. 
However, the inner Chandra ring does not show such obvious 
asymmetry, even though it appears to be spatially colocated 
with wisp 1. If this apparent symmetry is real, our identifica- 
tion of the first ion turning point with the Chandra ring ob- 
viously conflicts with our interpretation of the optical wisps. 
We believe that this discrepancy is another manifestation of 
the one-dimensional model's limitations. If the Chandra ring 
forms in the equatorial current layer, which is embedded in 
a higher latitude (|A| > 10°) flow (Begebnan 1999), perhaps 
having higher a (Arons 2004), the flow right in the equator 
might have lower velocity than the higher latitude compo- 
nent, with the optical wisps being features stimulated in the 
higher latitude flow by the large Larmor radius ions. Such a 
multi-dimensional environment, which is impUed by the ions 
being a net electric current, then might account for Doppler 
boosting of the optical wisps along with smaller or negligible 
Doppler boosting of the X-ray features. Testing dynamical 
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viability of such a scheme awaits the development of a two- 
dimensional, axisymmetric model, whose MHD aspects are 
likely to be similar to the nebular model recently proposed by 
Komissarov and Lyubarsky (2003). 

5.3. Nonthermal Particle Acceleration and Radiation 

Although our dynamical model and synthetic images as- 
sume a thermal (essentially, a monoenergetic) pair distribu- 
tion in the radiating plasma, they can be extended to the 
regime where the distribution is nonthermal, as is the case 
for the X-ray and optically emitting particles. The com- 
pressions in the field and plasma density associated with the 
shock structure should cause higher synchrotron emission for 
all electrons and positrons with Larmor radii smaller than 
the characteristic wisp separation, or for 7-1- < lO'*^, E± < 
5x10'^ eV, emitting at 160 keV < e„„ < 50- 100 MeV in the 
magnetic field whose strength increases from ^ y/o'Bnebuia ~ 
16 fiG at the leading edge of the shock structure to ^ 300 /xG 
at and beyond the outer reaches of the X-ray torus (Kennel and 
Coroniti 1984a; de Jager and Harding 1992). Thus, the wave 
motions seen in the Chandra band should display behavior 
similar to the optical wisps; that is, the features emitted from 
the X-ray inner ring discovered by Weisskopf et al. (2000) 
display similar variability as the optical moving wisps, and 
all propagate toward and through the X-ray torus. The obser- 
vations do show such resemblances, at the qualitative level. 

We have identified the Chandra X-ray ring as the site of 
the limit cycle in the unstable, reflected ion stream, not as 
the leading edge reverse shock in the pairs. The main rea- 
son to do this is the morphological similarities between the 
observed ring and our synthesized image. Whether the ac- 
tual pair shock is located at the Chandra ring as concluded 
by Hester et al. (2002) or interior to it depends on the shock 
acceleration mechanism for the pairs. Below we consider the 
consequences of several shock acceleration schemes. 

The pair shock heats the pair plasma to characteristic par- 
ticle energies 7wmd»J±c^ « 10'^7„6.5 eV, where 7^6.5 = 
7iwnrf/10^'^. The downstream magnetic field at and imme- 
diately behind the pair shock is weak. For rr « 3 x 10"^', 
B « \6Bs\(i{r/Rs±) /xG. Pairs compressed and heated by the 
pair shock radiate in the optical and infrared, with charac- 
teristic frequency I'c = 6 x lQi^'^^l^{r/Rs±) Hz. However, the 
synchrotron radiation time of newly shocked pairs at and just 
behind the pair shock is long, T,± w (3800/b2^^76.5)(/^.v±A)^ 
yr, far larger than the local flow time Taow = r/v{r) = 
0.5(/?s±/10^^^cm)(r//?j±)^ yr. Thus simple shock heating 
of the pairs implies that optical emission does not start un- 
til r > 10"* cm from the pulsar, which is larger than the radius 
of optical emission onset. 

Gallant et al. (1992) demonstrated that the pair shock can 
indeed thermaUze the upstream flow. The pair shock may, of 
course, also be a nonthermal particle accelerator, as has been 
assumed in all MHD models for the excitation of the Crab, 
starting with Rees and Gunn (1974). Diffusive Fermi acceler- 
ation (DFA) in transverse relativistic shocks is often invoked 
as the mechanism for such acceleration (see Gallant 2002 for 
a review). Since the energy gains per encounter of a particle 
with the shock are AS ^ E, this process might turn the pair 
shock into a nonthermal X-ray emitter, if the cross field diffu- 
sion time is short in the downstream medium, comparable to 
the cyclotron time ("Bohm diffusion") - such diffusion must 
occur, if the requisite multiple encounters are to exist. Diffu- 
sion speeds exceeding the downstream flow speed (~ c/3) re- 
quire quite strong downstream turbulence, however (Bednarz 



and Ostrowski 1996), probably at levels higher than observed 
in the two-dimensional pair shock simulations of Gallant et al. 
(1992). In the absence of three-dimensional simulations that 
allow evaluation of the turbulence level and that also aUow for 
cross field diffusion, DFA at the pair shock is not excluded 
as the means of creating nonthermal emission, although the 
constraints on the required diffusion rates clearly are rather 
demanding. 

Progressive acceleration with increasing radius, due to the 
resonant cyclotron absorption by the pairs of the high har- 
monic ion cyclotron waves emitted by the ions (Hoshino et 
al. 1992; E. Amato & J. Arons, in preparation), is an alter- 
nate to the commonly assumed DFA that has a natural setting 
in the context of the ion-doped relativislic shock whose spa- 
tial structure has been modeled here. The acceleration time 
is approximately the ion cyclotron time, which leads to en- 
ergy gains such that X-ray emission from the Chandra ring 
requires placing the now unobserved pair shock at a radius 
between 2 and 3 times smaller than the radius of the X-ray 
ring. A fuU discussion of this distributed acceleration model 
will be given elsewhere. 

5.4. The Acceleration of the High Energy Ions in the Wind 

So far we have not addressed the issue of how the ions are 
energized to the inferred values of 7 ~ 10*. One possibility is 
that the high energy of the ions and the high 4-velocity of the 
pulsar wind are due to the gradual "surf-riding" acceleration 
in the inner wind of the pulsar. If the wind begins near the 
Ught cylinder with ctlo = 1lo{Rl) ^ 1. where 7z, = 7(/?z,) is 
the Lorentz factor of a particle injected into the wind at r = R^, 
then particles riding on the field lines of the wind move on 
almost radial paths and increase their energies linearly with 
radius (Buckley 1977). This acceleration continues until the 
almost force-free, multidimensional flow in the inner region 

1/3 

ends at the radius Rff, which is not greater than Rl, where 

cr(r) drops to and 7wmd ~ o]^ (Beskin, Kuznetsova, & 
Rafikov 1998; J. Arons 2004, in preparation). Here Rjj is 
the radius where the flow speed first reaches the velocity of 
the fast magnetosonic wave; beyond this radius, the magnetic 
stress no longer can overcome the longitudinal inertia P7^,^ 
of the plasma, where p is the rest mass density. The acceler- 
ation in the iimer region has its origin in the particles being 
stuck to the field lines, with the balance of electromagnetic 
forces causing the field lines to move with an E x B velocity 
that is slightly below c, but increasing toward c with distance. 
This effect results in jmc^ increasing Unearly with radius for 
all particles in the wind, when r < Rff. For r > Rff, 'jwind 
stays approximately constant in laminar, ideal MHD relativis- 
tic flow. This conclusion has been formally demonstrated only 
when the poloidal field is monopolar, but is likely to be a ro- 
bust conclusion in practical terms from all relevant magnetic 
geometries (Begelman and Li 1994; J. Arons 2004, in prepa- 
ration).' 

' Early one-dimensional MHD models (Michel 1969; Goldreich & Julian 

I /3 

1970) found f cr^jj' as r ^ oo, a result that is an artifact of their one- 
dimensional (in effect, cylindrical) flow assumptions. The fact that Rff < oo 
in the multidimensional model of Beskin et al. (1998) results from their incor- 
porating the transverse (to the radial flow direction) electromagnetic stresses, 

while their result that Rff = (t^J^Rl 2> Rl is an artifact of using a strict 
monopole for the poloidal magnetic field. More physical models that have 
modest departures of the poloidal field from the monopole have Rff at most 
a few times Rl- Contopoulos and Kazanas (2002) are incorrect in their claim 
that the force free region extends to ~ (TioRl with 7 — » auo - the accelera- 
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Such a linear accelerator in the inner wind, an idea recently 
resurrected by Contopoulos and Kazanas (2002), reaches 
finite-inertia saturation too soon to provide sufficient accel- 
eration to the wind and the ions in the Crab. Our results for 
the particle content of the wind's equatorial sector with lati- 
tude range |A| < 10° suggest a Goldreich-Julian flux of heavy 
ions Ni = Ncj/Z = IQ'^^i/Zec w 2 x lO^V-^ s"', and a flux 
of electron-positron pairs N++N- = 2k±Ngj ~ lO'^Ni ^ 10^^-^ 
s"' (Harding & Muslimov 1998; Hibschman & Arons 2001; 
Mushmov & Hai-ding 2003). Then 



ctlo ^ 



Bl 



1 Ze^ mi 



(7) 



47r[m,w,(7?L) + 2ot±«±(/?l)]c^ 2 ot,c^ '«(>// ' 

Meff = Anip + 2K±Zm± . (8) 

For Crab pulsar parameters (magnetic moment fi ~ 10^°'^ cgs, 
n Ri 200 s"' , $ w 4 X 10'* V) and equatorial plasma properties 
inferred from our model of the wisp and ring dynamics (Z = 1 
or 2, nii = Anip with A = 1 or 4 and 2k± ~ 10"*), 

cTio « (1.3-1.7) X 10" > 1. (9) 
Our model for the wisps' dynamics suggests 7,0,, ^ 

0. lZe$/m,c^ au) ^ ct^q^. Since laminar acceleration of the 
whole relativistic wind does not explain the inferred 7,0,, w 
Jwind ^ 10" or a^indirs) ^ 10"^'^ in the equatorial flow (the 
only latitudes where we have substantial quantitative esti- 
mates for 7 and a in the wind), some additional physics must 
be at work to accelerate and demagnetize the wind in the equa- 
torial sector. 

The fact that the wisps (and the X-ray torus) lie in an angu- 
lar sector near the pulsar's rotational equator and that these re- 
gions appear to accumulate the majority of the instantaneous 
rotational energy loss offers an important clue to the addi- 
tional acceleration physics, a clue replicated in an increas- 
ing number of other young PWNs. That clue suggests that 
the equatorial region of this relativistic outflow has special 
properties not shared with the higher latitude regions, perhaps 
arising from the electric (return) current flowing in a crinkled 
sheet structure in an equatorial angular sector Substantial the- 
oretical evidence has accumulated to suggest that the ct ^ 1 
regions of pulsar winds have structure close to that of a split 
monopole for r larger than a few times Ri. This result is 
quite certain for the aligned rotator, both from the asymptotic 
analysis (Michel 1974) and from the self-consistent numeri- 
cal solution recently obtained by Contopoulos, Kazanas and 
Fendt (1999). While there is no similar demonstration that the 
asymptotic structure of the force-free oblique rotator's wind 
must be close to that of the oblique split monopole, whose 
field structure was recently obtained by Bogovalov (1999), 
we regard Bogovalov's model to be an excellent representa- 
tion of the likely structure of the electromagnetic fields in the 
wind of an oblique rotator outside the light cylinder and will 
assume those fields in this discussion of wind acceleration. 

In the equatorial sector < A < /, with / the inclination an- 
gle between the magnetic moment and angular velocity vec- 
tors, the electromagnetic fields have the form 



M 

Br=-f 



-{(t>-nt) 



M&mO 
E0=Bjy, Bg = 



r 

Wl~ 

Er = E, 



{(t>-VLt) 
, = 0. 
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tion dynamics is well described by force-free theoiy in a much more limited 
region than they assumed. 



Fig. 4. — Frozen-in cun'ent sheet structure of the oblique split monopole 
in the inner wind. Left: Meridional cross-section of the poloidal field struc- 
ture, showing the crinkled current sheet. The thick arrow indicates the in- 
stantaneous direction of the magnetic axis. Rotation is around the vertical 
axis. Right: Intersection of the current sheet with the equatorial plane. The 
toroidal magnetic field forms stripes with opposite directions between each 
current layer, as indicated by the aiTows between the cun'ent sheets. 



Here r, 6*, are the spherical radius, colatitude (measured from 
a z-axis with positive z along the neutron star's angular ve- 
locity vector), and rotational azimuth, respectively, measured 
from an arbitrary x-axis with increasing in the direction 
of rotation. Furthermore, / = ± 1 is a step function of pat- 
tern phase (r/pRi) + (t)-ilt, where c/3 is the radial outflow 
velocity. The fields reverse in sign every wavelength [at 
r = (n+ l)7r/3/?/^/2], while the magnitudes of the fields remain 
the same as those of the aligned split monopole. The effective 
monopole moment is M = kfji/Ri, with k a constant on the or- 
der of unity; A; = 1 .36 in the Contopoulos et al. (1999) solution 
for the aligned rotator A crinkled current sheet separates the 
regions of oppositely directed field, shown in Figure |4] (fol- 
lowing Bogovalov 1999). As is the case with the Deutsch 
(1955) fields of the vacuum rotator, these fields form a spi- 
rally wound pattern that corotates with the star Outside the 
equatorial sector, / = constant. For / < 7r/2, as is thought to 
be the case for the Crab pulsar (Yadigaroglu & Romani 1995), 
/ = 1 for < (7r/2) - / and / = -1 for (7r/2) + / < < tt. 

In the global, ideal model of the current flow, the ions re- 
side in the crinkled current sheet, which is frozen into the flow. 
Therefore, they participate in the same acceleration to which 
the plasma outside the current sheet is exposed. In the global 
wind models the current sheet has infinitesimal thickness. 
However, in reality the current flows in a crinkled layer of fi- 
nite thickness and may be subject to dissipative processes that 
can lead to the release of magnetic energy in the equatorial 
sector in the many decades in radius between r = Rff < IQ^Rl 
and r = = IO'^Rl. The dissipation of magnetic energy can 
lead to acceleration (and possibly heating) of the plasma (in- 
cluding the ions) to higher 4-velocities than is possible in the 
ideal flows so far studied. Whether such dissipation is induc- 
tive destruction of the current sheet, modeled as a sheet pinch 
(Coroniti 1990; Michel 1994; Lyubai'sky & Kirk 2001, Kirk 
& Skjaerassen 2003), or has something to do with the dissipa- 
tion of relativistically strong electromagnetic waves (Melatos 
1998), perhaps generated by instability of the current sheet 
and accelerating the ions through ponderomotive wave pres- 
sure (Arons 2003, 2004), remains an open question. 

6. CONCLUSIONS 

We have studied the internal structure and time variability 
of a collisionless shock in an electron-positron pair plasma 
with an energetically significant admixture of ions. Upon 
crossing the shock the ion component undergoes relativistic 
cyclotron instability that renders the shock structure very dy- 
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namic. The first loop of the E x B ion orbit achieves a limit- 
cycle behavior and acts as a quasi-periodic wave emitter with 
a period of half a local Larmor time of the ions. This time- 
dependence led us to apply our results to the collisionless ter- 
mination shock in the inner Crab Nebula where we identify 
the wisps with turning points in the drift orbit of the ions. The 
limit cycle of the instability produces large amplitude magne- 
tosonic waves on a time-scale of several months that proceed 
to move downstream from the shock. This time-scale is intrin- 
sic to the model and associates the appearance of the wisps to 
the mechanism of variability within the shock structure. This 
is in contrast to other models that generally do not explain the 
origin of the wisp perturbations, only their subsequent evo- 
lution. Our model naturally makes the region between wisp 
1 and wisp 2 in the optical (or the iimer Chandra X-ray ring 
and the torus) behave like a near zone of a radiating antenna 
and prescribes the pattern of brightness fluctuations to wisps 
1 and 2 and the moving wisps. This pattern has both a station- 
ary component oscillating in brightness (inner Chandra ring 
and the torus region) and moving wisps that cross the region. 
That is why the snapshots of the inner nebula taken with in- 
sufficient temporal sampling register motion yet on average 
always see two major wisps that do not seem to leave. 

The morphological agreement between the behavior of the 
wisp region in the Crab and the results of our simulations is 
an important argument in favor of the presence of ions in the 
pulsar outflows, which are often presumed to consist only of 
electron-positron pairs. Such ions are a natural candidate for 
the return current flow of the pulsar that starts at the auroral 
boundary of the polar cap and is then mapped into the equato- 
rial current sheet. In order to get a fit to the nebular dynamics 
we require an ion injection rate of IQ-''*^ s"'. The associated 
current is of the same order of magnitude as the Goldreich- 
Julian current of the Crab pulsar. Thus, by studying the in- 
teraction of the pulsar wind with the nebula we are able for 
the first time to get a window on the enigmatic electrodynam- 
ics of pulsars. Electrostatically, ions are the preferred return 



current carrier for pulsars with fl - fi > 0, which is thought to 
be the geometry applicable to the Crab pulsar (Yadigaroglu & 
Romani 1995). This means that for a rotator with antiahgned 
rotation and magnetic axes, the equatorial current would be 
carried by electrons, and their effect on the nebula might be 
different from the Crab case. 

The model of the time variability of collisionless shocks 
presented here is fairly general and can be applied to other 
PWNs, as well as other astrophysical sources where colh- 
sionless shocks arise. For the pulsar B 1509-58 and the as- 
sociated nebula Gaensler et al. (2002) estimate a variability 
period of on the order of 3-5 yr For Vela the spacing of 
semicircular X-ray arcs implies a period of variability sim- 
ilar to the Crab, yet there are no similar wisps in that rem- 
nant. This could be explained if Vela rings are not features 
in the equatorial flow, but rather in the polar direction (sim- 
ilar to the Crab halo (Hester et al., 1995)). The differences 
might also be a consequence of the later evolutionary state 
of the Vela PWN, where the environment has a chance to in- 
fluence the ternaination of the pulsar's wind (Chevalier 1998, 
2004). This underscores how much we are dependent on a re- 
liable deconvolution of the geometry and the evolutionary his- 
tory in all attempts to characterize the physics of non-thermal 
nebula energization by central compact objects. For sources 
other than the PWNs, colUsionless shocks are common in rel- 
ativistic AGN jets and thought to occur in gamma-ray bursts. 
Both settings are known to be quite variable, and whether ion- 
cyclotron instability can operate in these shocks remains to be 
studied. 
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APPENDIX 

APPENDIX: MODEL DERIVATION AND SIMULATION DETAILS 

Evolution equations for a time-dependent shock in the pair-ion system can be obtained from the divergence of the full energy- 
momentum tensor T"^ = nimiC^ufuf+wu^u^+P^^ + T^^. Here nitriiC^ is the rest energy density for ions, u = (j,'jf3) and 
Ui = (7( ,7!/3i) are 4- velocities of the pair and ion components, respectively, w is the enthalpy of the pairs, P"^ is the pressure 
tensor, and 7^ is the electromagnetic stress-energy tensor We assume that ion fluid is cold compared to pairs; thus, only 
pair pressure enters into f"'^. Beyond the shock the dynamics of the pairs is initially constrained to be in the plane orthogonal 
to the toroidal magnetic field 5^; therefore, we should allow for potentially anisotropic pressure: P" = -P±, P" = P^^ = P±, 
p(p(p - Allowing only for spatial variation in r for all quantities, the continuity equations together with the momentum and 
energy conservation laws and induction equation in spherical coordinates around the equatorial plane are then: 



4-M±l) + \-^{r'n±jpr) = 0, (Al) 
Oct r^ or 

d 1 ^ 

^Ma) + -jrir\ll3ir) = 0, (A2) 
Oct r^ or 

■S-inniMiC^Ui, + W--IU, + -^EgB^) +\-^ [r^iiiimiC^iuirf + w{urf - ^ {E^ 
Oct Ait r^ or 87r 



1 d P±-P\\ 1 , 9 

-r —PjL + - -inimiiuier+wiugY) - —! 

or r r 47i 

d , i „ „ . ^ d ^ 0, 1 



-E^e -Bl))] + - —P^ + ^ - -(nmiiutef+wiuef) - = 0, (A3) 

^ r"^ or r r Anr 



-(ni"/imiC^Ui0 + w^U0 - —ErB^) + -5- tt [r^inimic\rUi0 + wUrUg - -^ErEg)] 
Oct A-K r^ or 4tt 

1 2 1 

+ -[nimiC UirUiO + WUrUg - —ErEg] = 0, (A4) 

r 4w 
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d \ \ d 1 

— [niniiC^-iJ + W7^ -P^+—{E^ + B^)] + ^Tr [r^ininiiC^riUi, + w-fUr + —EgB^)] = 0, (A5) 
Oct Svr Or An 

d 19, BAf3r 

—B^ + -—{r'B^P,) ^ = 0. (A6) 

Oct r'^ or r 

In jA6> we used the MHD condition E + (3 x B = 0. It can be shown that for variations on the ion cyclotron time-scale and 

moderate harmonics thereof, departures from ideal MHD in the pairs are negligible. 

We proceed to normalize all dimensional quantities in these equations in terms of the upstream ion gyration radius and cyclotron 

time as in ^ Using the normalized quantities and the generalized adiabatic index for anisotropic pressure as defined in ^ the 

enthalpy can be written as w = p(l + r/(r- l)P±/p) = ph , where p is the rest energy density of the pairs, and h is the specific 

enthalpy. Using the equations of motion for the ions, 

^-i^ = ?i[i?+^xi3], (A7) 
act r 7; m,- 7,- 

duis luigUir Ze Uj 

-J— + = —{E+—xB]g, (A8) 

act r 7; nii 7/ 

we can rewrite the normalized pair equations ( IA1IA6> as follows: 

o \ c) 

j-{pi)+-j-(r''p-il3r) = Q, (A9) 
at or 

^(phj^P,. + a±prBl)+^j-y(phj^Pj + ^a±(l-pj + P^)Bl)] 

+4 ^Pi- + ^^^-^ - -phl^Pl - ^± -PlBl + eNiiPie - Pe)B^ = 0, (AlO) 
r^ or r r r ^ 

^(phj^lSe + a±PeBl) + [,^(ph-f^f3rPg + a±/3,/3eB^)] 
+^ [phj^Me + cT±MeBl] + eNi(A - l3ir)B^ = 0, (Al 1) 

|- [phj^ -P±+\<J±il + /3^ + pj)Bl] 

1 _9 

■ 9r '" 

d 19, Ba 
-B,.-,-yB,p,.)-^^ 

where energy fractions a±, cr, and e are defined in (|5Jl. For flows with low magnetization such that a± 0, ct, 0, but e 
finite, these equations reduce to the system (I1I4> . The Poisson equation and Ampere's law can be written in terms of normalized 
quantities as 

p± + Pi = <TiV-E, (A14) 

j±+j,. = a,(VxS-|-S), (A15) 

where the charge densities of the pair fluid (p±) and the ions {pi) are normalized in terms of the upstream ion density A^i,, and 
similarly for the currents j^, jj. When cr, is small, field gradients are on scales commensurate with r,i in all directions, and time 
evolution is on scales lo~Iy, the expressions on the right-hand side of ( IA14I IaTsI can be neglected and the flow is charge and 
current quasi-neutral.^ 

We have developed a code that solves the full equations (IA9IA13t : however, for the region inside of the X-ray torus in the Crab 
Nebula, which is the region of interest for the wisp dynamics, the inclusion of the magnetization does not change our results. It 
is important to include these terms to study the evolution of the flow in the outer nebula where the magnetic field is compressed 
to equipartition, so in the f ollowing we describe the structure of the full simulation. 

We discretize eqs. (I1I4> on an equidistant grid in radius ranging from 0.1 to lOrm with 1000 cells. Our simulation domain 
contains both the regions of supersonic flow upstream of the shock and the subsonic downstream region, while the position of the 
shock is determined dynamically. To advance the pair equations, we use the Lax-Friedrichs central-difference forward-in-time 
scheme. Despite its simplicity, it handles well both the upstream and the downstream regions, including wave propagation and 
scattering off the shock. To maintain numerical stability of a centered scheme, the method introduces inherent diffusivity, which 
spreads the reverse shock over 10 cells. This is not a concern since the shock is still very well localized on the grid, and the 
wave amplitude diffusion is less than 10% over the simulation domain. At every timestep we advance the following conservative 
variables: D = 7/), F = phj-f3r + o'±f3rB^, G = phj^(3g + a±PgB^, and r = phj--P± + \(t±{1+ f3^ + (3j)B^. These variables need 
to be converted to the primitive variables p, P±, jSr.g, and B^. We use a numerical root-finding routine to solve for pressure P± 
from P± = (r- l)e,, where the energy density is given by 

e* = [r + D(l-7,) + Pi(l-72)-a±B2(2-l)]l. (A16) 

7* li 

^ We are indebted to Z.-Y. Li, who first pointed out this possible simplification to us. 



+^ — [r\ph-f^f]r + a±prBl)] + emfiief3r-PirPe)B^ = 0, (A12) 
r-^ Or ^ 

1 d , BsB,- 

-B^ + -—(r^B^f3,) ^ = 0. (A13) 
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At every iteration of the root finding the Lorentz factor 7* is expressed as a function of r, F, G, P±, and B0 through 



(r + Pj. + ia±B2/72)2-(F2 + G2)' 



which can be reduced to a cubic in 7^. Having found P± and 7, pair flow velocities are obtained from 

Pr= ^, (A18) 

Pe = ^. (A19) 

r+P± + {<J±% 

The initial state for the calculation is the steady state strong MHD shock. We construct this state by first computing the jump 
conditions between the upstream (subscript 1) and downstream (subscript 2) side of the shock from: 



Pi7?/3h + :t<^±(1 + l3n)Bli = PihiljP^i + ^'^±(1 + P^2)Bl2 + P^2, (A21) 



PlllPrl = PlllPrl, (A20) 

ia±( 1 + l3\)Bl, = P2h2llf3j2 + 

PX-ll(3r\+(7±fir\Bli=P2h2llPrl + (y±Pr2Bl2. {All) 
B,p\(3r\=B^2f3r2, (A23) 

where the normalized quantities on the upstream side are given by 71 = 10^, B^i = I, pi = the flow is cold, P±i = 0, and is 

assumed to have no ^-component of velocity. Although there exist analytical approximations to the jump conditions for various 
limiting values of a± (Kennel and Coroniti 1984a), we solve the general case numerically. Having found the jump conditions, 
we select the radius of the shock rj (typically, = Irm) and solve for the steady state flow from conservation laws ( IA9IA12t 
with the time derivatives omitted. 

The flow solution can be expressed in terms of constants ci = r^jpPr from the equation of continuity, 02 = r^phj^(3^+r^a±PrB^ 
from the energy equation and C3 = rB^fi,- from the induction equation. These constants can be evaluated at the shock and then 
used together with the equation of state to express P^, p, and as a function of /3r. When substituted into the momentum 
equation 



l^^[r\ph^^(3l + \a±{\- (3l + 0l)Bl)] + ^^P^ + ^^^-!^ =Q, (A24) 
or 1 ^ n or 



we get a differential equation for (ir, which is integrated numerically inward and outward from the shock to obtain the flow 
solution. 

The dynamical evolution equations JA9IA12t are driven by the currents induced in the pairs by the motion of the ions. To 
compute the ion quantities A^, and /3,>,e, we use the 1 .5 -dimens ional particle-in-cell method. Ions are represented as macroparticles 
whose velocities are calculated by advancing eqs. ( IA7IA8> . The fields are extrapolated from the grid to the positions of the 
particles, and after advance the particle densities and currents are deposited at gridpoints. In both of these operations we utilize 
volume weighting to properly account for the diverging spherical geometry. 

The simulation domain is augmented by a buffer zone of width equal to 1 ru\ where the field is kept constant. This zone is 
needed to accurately represent the contribution of ions that would normally leave the domain and then reenter the outer edge 
of the grid as part of their E x B trajectory. The outer boundary conditions are designed to mimic the rest of the nebula that 
is not simulated and to reduce the unwanted reflections. Two methods are used for this. For simulations with / = ruxjrs <C 1, 
or essentially planar geometry, the reflections off the right edge in subsonic outflow can be minimized by applying the method 
of characteristic decomposition (Thompson 1987). The method works by explicitly setting the left-going characteristics to zero 
at the boundary. For simulations with significant spherical geometry effects (/ ^ 1) less complicated boundary conditions of 
variable extrapolation produce acceptable results. 
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